function Rates = MISO_compute_rate(qbit,Nsub,L,TRX,fc,BW_des,Ptx,BS,showfig)

BSN1 = BS.n1;
BSN2 = BS.n2;
BSD1 = BS.d1;
BSD2 = BS.d2;







Ngrid = 1000;
% L = 15e-1;
% 
% factor = 0.3;
% x1 = 40e-1*factor;
% y1 = 0e-1*factor;
% z1 = 90e-1*factor;
% x2 = -18e-1*factor;
% y2 = 0e-1*factor;
% z2 = 20e-1*factor;
x1 = TRX(1);y1 = TRX(2);z1 = TRX(3);
x2 = TRX(4);y2 = TRX(5);z2 = TRX(6);
z12 = z1^2;z22 = z2^2;
% 
% fc = 30e9;
c = 3e8;
% BW_des = 1.5e9;
% Ptx = 0;%dBm


Ptx_lin = 10.^((Ptx-30)/10); % in watts
Lb = c/fc;
spac = Lb/2;
hp = 6.625e-34;
kb = 1.38e-23;
Tk = 290;


LIS_f = 0;
N = round(L/spac);


% RIS configuration
xcoord=linspace(-L/2,L/2,N);
ycoord=xcoord;
[xx,yy]=meshgrid(xcoord,ycoord);
coord_mat=xx+1i*yy;
a = (sqrt((xx-x1).^2+(yy-y1).^2+z1^2)+sqrt((xx-x2).^2+(yy-y2).^2+z2^2))/2;
a_vec = linspace(min(a(:)),max(a(:)),Ngrid)';% column vector
Distmat = a*2;
distgrid = a_vec*2;
Distmat1 = sqrt((xx-x1).^2+(yy-y1).^2+z1^2);
Distmat2 = sqrt((xx-x2).^2+(yy-y2).^2+z2^2);

BScoord1 = spac*(-BSN1/2+1/2:BSN1/2);
BScoord2 = spac*(-BSN2/2+1/2:BSN2/2);
[c1,c2] = meshgrid(BScoord1,BScoord2);
BSx = c1*BSD1(1)+c2*BSD2(1)+x1;
BSy = c1*BSD1(2)+c2*BSD2(2)+y1;
BSz = c1*BSD1(3)+c2*BSD2(3)+z1;

[ModBFmat,QBFmat,modfn] = Beamforming(qbit,L,TRX,fc,BW_des,Ptx,showfig);

hcf = zeros(numel(c1),Nsub); % cascaded channel h_c(f)
WBS_PS = zeros(numel(c1),1);% weighting matrix for BS
WBS_TTD = zeros(numel(c1),Nsub); % frequency-independent BF
freqvec=2*pi*linspace(fc-BW_des/2,fc+BW_des/2,Nsub);

% gpuidx = getgpu(8e9);
gpuidx = 0;
% disp('gpuidx:')
% disp(gpuidx)
if gpuidx

    gpuDevice(gpuidx);
    xx =    gpuArray(complex(xx));
    yy =    gpuArray(complex(yy));
    % zz =    gpuArray(complex(zz));
    x2 =    gpuArray(complex(x2));
    y2 =    gpuArray(complex(y2));
    z2 =    gpuArray(complex(z2));
    BSx =    gpuArray(complex(BSx));
    BSy =    gpuArray(complex(BSy));
    BSz =    gpuArray(complex(BSz));
    freqvec  = gpuArray(complex(freqvec));
    c = gpuArray(complex(c));
    hcf = gpuArray(hcf);
    WBS_PS = gpuArray(WBS_PS);
    WBS_TTD = gpuArray(WBS_TTD);
end

for nn = 1:numel(c1)

    Distmat = (sqrt((xx-BSx(nn)).^2+(yy-BSy(nn)).^2+BSz(nn)^2)+sqrt((xx-x2).^2+(yy-y2).^2+z2^2));
    Distmat1 = sqrt((xx-BSx(nn)).^2+(yy-BSy(nn)).^2+BSz(nn)^2);
    Distmat2 = sqrt((xx-x2).^2+(yy-y2).^2+z2^2);
    antt_delays=Distmat(:)/c;
    AF_product=antt_delays*freqvec;
    h_0 = (Distmat1(:)*freqvec/c).*(Distmat2(:)*freqvec/c);
    Wb_mch=exp(-1i*AF_product)./h_0;
    hcf(nn,:) = ModBFmat(:).'*Wb_mch;
    WBS_PS(nn) = sum(hcf(nn,:));
    WBS_PS(nn) = conj(WBS_PS(nn))./abs(WBS_PS(nn))/sqrt(numel(c1));
    WBS_TTD(nn,:) = conj(hcf(nn,:))./abs(hcf(nn,:))/sqrt(numel(c1));

end

  SCpowers_PS = abs(WBS_PS.'*hcf).^2;
  SCpowers_TTD = abs(sum(WBS_TTD.*hcf,1)).^2;

% ideal: N basestation antenna with 0 spacing. ideal gain.
    Distmat = (sqrt((xx-x1).^2+(yy-y1).^2+z1^2)+sqrt((xx-x2).^2+(yy-y2).^2+z2^2));
    Distmat1 = sqrt((xx-x1).^2+(yy-y1).^2+z1^2);
    Distmat2 = sqrt((xx-x2).^2+(yy-y2).^2+z2^2);
    antt_delays=Distmat(:)/c;
    AF_product=antt_delays*freqvec;
    h_0 = (Distmat1(:)*freqvec/c).*(Distmat2(:)*freqvec/c);
    Wb_mch=exp(-1i*AF_product)./h_0;
    SCpowers_0spacing = numel(c1)*abs(ModBFmat(:).'*Wb_mch).^2;

    SCpowers_ideal = ones(size(freqvec));
    Pow = sum(modfn.^2)*spac^-4*(a_vec(end)-a_vec(1))/Ngrid*c/(freqvec(2)-freqvec(1))*pi./h_0(end/2,end/2)^2;
    SCpowers_ideal = numel(c1)*SCpowers_ideal*Pow/sum(SCpowers_ideal);

  
        %% Compute rate
        
        Delta_f=freqvec(2)-freqvec(1);
        phin_nr=hp*freqvec; % ε=hν
        phin_dr=exp(hp*freqvec/(kb*Tk))-1;% Bose distribution factor
        
        phin=phin_nr./phin_dr;  % in watts
        
Rates = zeros(3,1);

        Hwb_sq=SCpowers_PS;%
        N_eff=phin./Hwb_sq;
        mu_df=wfill(N_eff,Ptx_lin/Delta_f);
        powalloc=mu_df-N_eff;
        powalloc(powalloc<0)=0;
        SNR_eff=powalloc./N_eff;
        Rates(1)=Delta_f*sum(log2(1+SNR_eff));
        
        Hwb_sq=SCpowers_TTD;
        N_eff=phin./Hwb_sq;
        mu_df=wfill(N_eff,Ptx_lin/Delta_f);
        powalloc=mu_df-N_eff;
        powalloc(powalloc<0)=0;
        SNR_eff=powalloc./N_eff;
        Rates(2)=Delta_f*sum(log2(1+SNR_eff));
        % 
        Hwb_sq=SCpowers_0spacing;%
        N_eff=phin./Hwb_sq;
        mu_df=wfill(N_eff,Ptx_lin/Delta_f);
        powalloc=mu_df-N_eff;
        powalloc(powalloc<0)=0;
        SNR_eff=powalloc./N_eff;
        Rates(3)=Delta_f*sum(log2(1+SNR_eff));

        Hwb_sq=SCpowers_ideal;%
        N_eff=phin./Hwb_sq;
        mu_df=wfill(N_eff,Ptx_lin/Delta_f);
        powalloc=mu_df-N_eff;
        powalloc(powalloc<0)=0;
        SNR_eff=powalloc./N_eff;
        Rates(4)=Delta_f*sum(log2(1+SNR_eff));

        
        
        